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ABSTRACT 

The unusual flux variations of the pre-main-sequence binary star KH 1 5D have been attributed to occultations 
by a circumbinary disk. We test whether or not this theory is compatible with newly available data, including 
recent radial velocity measurements, CCD photometry over the past decade, and photographic photometry over 
the past 50 years. We find the model to be successful, after two refinements: a more realistic motion of the 
occulting feature, and a halo around each star that probably represents scattering by the disk. The occulting 
feature is exceptionally sharp-edged, raising the possibility that the dust in the disk has settled into a thin layer, 
and providing a tool for fine-scale mapping of the immediate environment of a T Tauri star. However, the 
window of opportunity is closing, as the currently visible star may be hidden at all orbital phases by as early as 
2008. 

Subject headings: stars: pre-main-sequence — stars: individual (KH 15D) — circumstellar matter — open clus- 
ters and associations: individual (NGC 2264) 
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1. INTRODUCTION 

After a newborn low-mass star emerges from the dust of 
its parent molecular cloud, it spends a few million years as 
a T Tauri star. This phase is characterized by optical vari- 
ability and chromospheric emission lines, as the star accretes 
gas, propels a bipolar outflow, and contracts onto the main 
sequence (see, e.g., the review by Bertout 1989). A conven- 
tional boundary line is drawn between "classical" and "weak- 
lined" T Tauri stars, depending on whether the equivalent 
width of Ha emission is greater or smaller than some thresh- 
old value, typically 5-10 Angstroms. Of course, this bound- 
ary line is somewhat arbitrary, and transitional cases are to be 
expected, but we are aware of only one T Tauri star, KH 15D, 
that seems to mock this definition by alternating between clas- 
sical and weak-lined with clocklike regularity (Hamilton et al. 
2003). 

This is not the only unusual property of KH 15D. Although 
it was classified as a variable member of the young cluster 
NGC 2264 by Badalian & Erastova (1970), the peculiarity 
of its variations was not appreciated until Kearns & Herbst 
(1998) found that the object fades by 3 magnitudes every 48 
days. Hamilton et al. (2001) and Herbst et al. (2002) pro- 
posed that the star was being periodically eclipsed by cir- 
cumstellar material. Further monitoring revealed qualitative 
changes in the eclipses, most notably a gradual increase in 
their depth and duration, as well as the progressive disappear- 
ance of a "re-brightening" event that once characterized the 
mid-eclipse phase. The dramatic growth in the equivalent 
width of Ha emission occurs during the transition from the 
bright state to the faint state. Apparently the occulting ma- 
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terial blocks the photosphere of the star and thereby reduces 
the contrast between the photosphere and the Ha-emitting re- 
gion: the star has a "natural coronagraph" (Hamilton et al. 
2003). Eclipse-related variations in both optical polarization 
(Agol et al. 2004) and molecular H2 emission (Deming, Char- 
bonneau, & Harrington 2004) have also been observed. 

For several years, the cause of the eclipses and of their evo- 
lution were the subject of much speculation (Grinin & Tam- 
bovtseva 2002, Herbst et al. 2002, Winn et al. 2003, Barge 
& Viton 2003, Klahr & Bodenheimer 2003, Agol et al. 2004, 
Barrado y Navascues et al. 2004). Recently it was proposed 
that KH 15D is an eccentric binary system that is being oc- 
culted by the edge of a circumbinary disk (Winn et al. 2004, 
Chiang & Murray-Clay 2004). In this scenario, we are view- 
ing the system nearly edge-on, and the sky projection of the 
disk acts as a dark screen that gradually covers the binary. 
The advance of the screen is a consequence of the nodal pre- 
cession of the disk, which is inclined relative to the binary 
plane. At present, one member of the binary is hidden at all 
orbital phases, and the other member undergoes an eclipse 
each time its orbital motion brings it behind the screen. Winn 
et al. (2004; hereafter, Paper I) showed that this model unified 
the diverse properties of KH 15D and provided a quantitative 
fit to the available measurements of the eclipses and their evo- 
lution. Chiang & Murray-Clay (2004) arrived at the same idea 
and argued that the disk must be warped and radially narrow 
(a "ring") in order to maintain rigid precession. 

It has since been confirmed that KH 15D is a spectroscopic 
binary with an eccentric orbit (Johnson et al. 2004). In addi- 
tion, a wealth of new photometric data has become available, 
from analyses of archival photographs (Johnson et al. 2005, 
Maffei, Ciprini, & Tosti 2005), and modern CCD observa- 
tions (Hamilton et al. 2005; Barsunova, Grinin, & Sergeev 
2005; Kusakabe et al. 2005). Our main motivation for the 
work described in this paper was to test whether or not the 
circumbinary-disk model is consistent with these new data. 
In §§ 2 and 3, we describe our compilation of the data, and in 
§ 4 we specify the model and its parameters. Section 5 is a di- 
rect extension of Paper I: we attempt to fit the observed times 
and durations of the eclipses, and the observed radial-velocity 
variations. 



2 



Winn et al. 



A second motivation for this work was to understand certain 
features of the data that had not previously been explained. 
What is the origin of the light that is received during eclipses? 
Why is the depth of the eclipses slowly increasing with time? 
What causes the nearly-repeatable flux variations that are seen 
throughout ingress, mid-eclipse, and egress? We found that 
by making a minor elaboration to the model — surrounding 
each of the stars with a faint blue halo — all of these obser- 
vations could be accurately described. This is shown in § 6. 

Finally, in § 7 we discuss the interpretation of the fitted 
model parameters, in particular those that describe the halos. 
The most likely interpretation is that scattering by the disk 
creates an apparent halo around each star. We suggest some 
future observations and calculations that may clarify the in- 
terpretation. If these are successful, and if the occulting edge 
is as sharp as we suspect, then it will be possible to use the 
occulting edge to map out the environment of the visible star 
with wondrously fine detail. This would broaden the signifi- 
cance of this system from being a fascinating puzzle to being 
a crucial fortuitous case that will provide a deeper understand- 
ing of all young stars. 

2. PHOTOMETRY 

Two types of photometry are available for this system: pho- 
tographic photometry, generally from years before 1985, and 
CCD photometry, from 1995 to the present. In all cases, the 
measurements refer to the total light from the system, because 
the system has not been spatially resolved (with the possi- 
ble exception of an H2 filament observed by Tokunaga et al. 
2004). The data are very heterogeneous, having been gathered 
with dozens of different telescopes and reduced with different 
procedures. This section describes the available data and our 
efforts to allow them to be modeled within a single frame- 
work. In general, our attitude was cautious, in the sense that 
we wished to maximize the reliability of the data even at the 
expense of discarding some of the data. 

2. 1 . Photographic Photometry 

As a fairly unobscured example of a young stellar cluster, 
rich with variable stars, NGC 2264 has long been a popular 
target for photometric campaigns. Many photographic plates 
of NGC 2264 taken over the past century have been preserved 
in observatory archives. Four analyses of archival observa- 
tions of KH 15D have been published: 

1 . Winn et al. (2003) used blue-sensitive exposures from 
the Harvard College Observatory to show that in the 
first half of the 20th century, the system was rarely 
(if ever) fainter than its modern out-of-eclipse level by 
more than 1 magnitude. Given the number of plates an- 
alyzed, they concluded that the system spent less than 
20% of the time in such a faint state. More accurate 
photometry was impossible because of the glare from a 
nearby star, HD 47887. 

2. Johnson & Winn (2004) performed profile photome- 
try on a digitized set of infrared-sensitive exposures 
from the 92/67 cm Schmidt telescope of the Astrophys- 
ical Observatory of Asiago, Italy, from the time period 
1967-1982. The glare from HD 47887 is less problem- 
atic on the infrared plates than on the blue plates. The 
relative magnitudes have a typical accuracy of 0. 1 mag. 
They were placed on the standard Cousins / band mag- 
nitude scale using a set of reference stars whose / band 



magnitudes had been measured with CCD observations 
by Flaccomio et al. (1999). The zero point uncertainty 
is 0.14 mag. 

3. Maffei et al. (2005) analyzed a collection of blue and 
infrared plates that were originally obtained by the 
lead author, P. Maffei, between 1955 and 1970 us- 
ing 5 different telescopes. Relative photometry was 
performed visually, with an estimated accuracy rang- 
ing from 0.05 mag to 0.2 mag depending on the tele- 
scope. The blue magnitudes were converted to the 
Johnson B band using a set of reference stars with 
photographic and photoelectric B magnitudes given by 
Walker (1956). The infrared magnitudes were con- 
verted to the Cousins / band, after adopting an /-band 
magnitude for each reference star based on an extrapo- 
lation of its B - V color. 

4. Johnson et al. (2005) gathered and digitized a collec- 
tion of 87 plates, representing 8 different telescopes and 
22 different filter/emulsion combinations, from the time 
period 1954-1997. The magnitudes of KH 15D were 
measured and reduced to UBVRI magnitudes using the 
procedure of Johnson & Winn (2004). 

The latter 3 samples have 56 plates in common, giving us 
the opportunity to check for consistency between the different 
photometric methods. The plates in common are all infrared 
plates from one of the two Asiago Schmidt telescopes: 38 are 
from the 92/67 cm telescope, and 18 are from the 50/40 cm 
telescope. Figure ^ shows a comparison of Im, the magni- 
tude reported by Maffei et al. (2005), and Ij, the magnitude 
reported by either Johnson & Winn (2004) or Johnson et al. 
(2005) based on the same plate. There are significant discrep- 
ancies. For the 92/67 cm plates, the relation between Ij and 
Im is well-fitted by a straight line, 

7,-13.80 = 0.67(4,-14.27), (1) 

as plotted in Fig. Q This relation may reflect the different 
methods for determining the zero point and for coping with 
the nonlinear response of the emulsions. For the 50/40 cm 
plates, there is no obvious relation between Ij and Im- The 
best-fitting line (the dotted line in Fig.[Q represents an anti- 
correlation between the two sets of results, which does not 
make sense, and suggests that the photometric uncertainties 
have been underestimated by one or both sets of authors. 
Of all the infrared plates, Johnson et al. (2005) found the 
50/40 cm plates to be the most difficult to analyze, because 
of the large plate scale and consequently poor sampling of the 
stellar images. 

In this work, we attempt to fit only the /-band data, mainly 
because the majority of modern CCD observations were made 
in the / band, and also because we judge the photographic 
photometry to be most reliable for the infrared plates. We 
briefly discuss the data at other wavelengths in § but we 
do not use those data to determine our model parameters. Be- 
cause we do not know how to resolve the discrepancy between 
the two different sets of /-band results, we do not consider the 
Asiago 50/40 cm data further. For the Asiago 92/67 cm data, 
we use Ij when it is available, and otherwise we place Im onto 
the Ij system 9 using Eq. Q. For the data from all of the other 

9 We transformed Im into Ij, rather than the other way around, because the 
zero point of the Ij scale was set with reference to CCD /-band observations, 
and we therefore expected Ij to provide better consistency with the modern 
photometry. 
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FIG. 1 . — Comparison between the /-band magnitudes of KH 15D reported 
by two different groups based on the same photographic plates. Each point 
represents a single infrared plate, with 1m defined as the value reported by 
Maffei et al. (2005), and Ij as the value reported by either Johnson & Winn 
(2004) or Johnson et al. (2005). Different symbols correspond to different 
telescopes. The solid line is a linear least-squares fit to the 92/67 cm data, 
and the dotted line is a linear least-squares fit to the 50/40 cm data. 



telescopes, we use the published results without modification. 
The resulting data set consists of 81 measurements. The data 
are plotted in Fig. |2] along with the CCD photometry from 
more recent times. 

2.2. CCD Photometry 

Over the last decade, KH 15D has been monitored inten- 
sively. The first few years of observations used the 0.6 m tele- 
scope at Van Vleck Observatory on the campus of Wesleyan 
University, in Connecticut, and were designed to monitor the 
variability of a large sample of stars in NGC 2264. After 
two seasons of observing, Kearns & Herbst (1998) identified 
KH 15D as an especially interesting variable star. Hamilton et 
al. (2001) focused specifically on KH 15D, presenting a third 
year of photometry along with optical spectroscopy and a dis- 
cussion of possible explanations for the peculiar eclipses. An 
observing campaign involving many observatories was orga- 
nized, and the first year's results were described by Herbst et 
al. (2002). Most recently, Hamilton et al. (2005) presented 
both new photometry and a comprehensive summary of 9 
years of CCD photometry involving a dozen different tele- 
scopes. We refer the reader to that work for a detailed discus- 
sion of the calibration and internal consistency of the various 
data sets involved. Most of the photometry is in the Cousins / 
band. A notable exception is the season of 2001-2002, when 
simultaneous Johnson VI measurements were made with the 
Yale 1 m telescope on nearly every clear night. There are also 
a smaller number of other observations in the UBVR bands. 

Barsunova et al. (2005) have independently monitored 
KH 15D in the Johnson V and Cousins R and / bands since 
2002, using the 0.7 m telescope of the Crimean Astrophysi- 
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cal Observatory. Eighteen of their / band measurements were 
taken within 6 hours of a measurement by Hamilton et al. 
(2005), which, as in the previous section, gives us an oppor- 
tunity to check for consistency between the teams' different 
photometric procedures. We found them to be consistent with 
the linear transformation 

I H = 0.97/ B + 0.27, (2) 

which is essentially a difference in the zero point of the mag- 
nitude scale. 

Recently, Kusakabe et al. (2005) presented near-infrared 
(JHK S ) photometry of KH 15D spanning two observing sea- 
sons, from 2003 to 2005. As stated previously, in what fol- 
lows we are concerned almost entirely with the /-band data, 
although we do remark on the near-infrared data in § Q 

For our model-fitting purposes, we considered all of the 
/-band data of Hamilton et al. (2005) and Barsunova et al. 
(2005), after transforming the latter data set according to 
Eq. 10 for the sake of uniformity. In addition, V. Grinin 
kindly provided us the results from the 2004-2005 observ- 
ing season that have not yet been published, which we also 
transformed according to Eq. (0 and added to the data set 
under consideration. In total, there are 6780 data points, of 
which 87 are from the Russian group. The uncertainties range 
widely from a few millimagnitudes, for some observations 
when KH 15D was in its bright state, to 0.25 mag, when it 
was in its faint state. 

2.3. Summary of the Photometry 

The entire /-band time series, from 1956 until the present, 
is shown in Fig. [2] Many interesting features of the system's 
photometric history can be seen in this figure. Before about 
1960, there is not much evidence for variability greater than 
0.1 mag. 10 Beginning in the early 1960s and lasting until 
at least the early 1980s, the system varied between / « 13.6 
and 14.5. By 1995, the variations were between 14.5 and a 
faint state that has gradually darkened from about 17 to 18.5. 
Although it is not apparent in Fig. [2] both the pre-1980 and 
post- 1995 variations are periodic with a 48.4 day period, and 
the phase of minimum light shifted by nearly 180 degrees be- 
tween those time periods (Johnson & Winn 2004). 

3. RADIAL VELOCITIES 

Possible time variations of the Doppler shift of starlight 
from KH 15D were first reported by Hamilton et al. (2001) 
and Hamilton et al. (2003). Subsequently, Johnson et al. 
(2004) showed that KH 15D is indeed a spectroscopic binary, 
based on a series of high-resolution optical spectra using tele- 
scopes at the Keck, Magellan, and McDonald observatories. 
With 16 measurements spanning 1.3 yr with accuracies rang- 
ing from 0.2 to 0.6 km s , Johnson et al. (2004) showed that 
the Doppler shifts could be explained as radial velocity varia- 
tions of a star in an eccentric Keplerian orbit. 

If the occulting edge of KH 15D is sharp enough to resolve 
the stellar photosphere, then the radial velocity measurements 

10 Maffei et al. (2005) found evidence for periodic variations in B between 
1958 and 1963, a time when /-band variations at the same level could be ex- 
cluded. Mo st of the relevant data is from the Asiago 50/40 cm Schmidt tele- 
scope. In ij |2.1l we showed that two different groups arrived at significantly 
different results when analyzing the same infrared plates from that telescope, 
suggesting that the uncertainties may be larger than they appear. The un- 
certainties should be still larger for the blue plates because of the difficulty 
presented by the neighboring bright star. Hence, although the possibility that 
the eclipses began in B prior to / is interesting and worthy of further investi- 
gation, we do not consider this possibility further in this paper. 
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FIG. 2. — Fifty years of Cousins /-ban d ph otometry of KH 15D. The post-1995 data have been averaged into 6 hour bins. (The dotted lines mark the time 
ranges that are plotted in the panels of Fig. 1121 



taken during ingress or egress are subject to systematic er- 
rors due to the Rossiter-McLaughlin effect (Rossiter 1924, 
McLaughlin 1924). This term refers to the apparent velocity 
perturbation of a partially occulted star due to stellar rotation. 
For example, if only the receding half of a star is visible, then 
the starlight is redshifted in excess of the shift one would ex- 
pect from only the star's center-of-mass motion. This effect 
has long been observed during eclipses of binary stars (see, 
e.g., Worek 1996, Rauch & Werner 2003), and more recently 
during transits of extrasolar planets (Bundy & Marcy 2000, 
Queloz et al. 2000, Snellen 2005, Winn et al. 2005). To esti- 
mate the magnitude of the effect, suppose that a star of radius 
R is occulted by a straight-edged semi-infinite screen, as de- 
picted in the left panel of Fig. [3] The minimum distance be- 
tween the edge and the stellar center is r. The angle between 
the normal to the edge and the sky projection of the stellar 
north pole is cf). The fraction of the total stellar flux (/) and 
the mean line-of-sight velocity (AV) of the exposed portion 
of the star are given by 



AV: 



4( 



COS u 



-uW- 



V rot sin h sine 



(l-u 2 ) 3 / 2 



-u\f\~- 



(3) 



(4) 



where u = r/R, V rot is the stellar rotation speed and 7* is the 
inclination of the stellar rotation axis with respect to the line 
of sight. These expressions neglect differential rotation and 
limb darkening. 

For the case of KH 15D, Hamilton et al. (2005) measured 
Vrotsin/* = 6.9 ±0.3 km s -1 for the currently visible star. We 
estimate <fi « 20°, using the model of Paper I, and assuming 
the stellar rotation axis to be perpendicular to the stellar or- 
bital plane. The right panel of Fig.[5]shows the resulting mag- 
nitude of AV as a function of /. For the perturbations to be 
smaller than 0.2 km s" 1 , and hence smaller than the measure- 
ment uncertainties in the radial velocity data, we need / > 0.9. 
Thus, for model-fitting purposes, we consider only the 12 ra- 
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FIG. 3. — The Rossiter-McLaughlin effect produced by a straight-edged, 
semi-infinite occultor. Left. — Diagram of the star, its projected rotation axis, 
and the occulting edge. Right. — The magnitude of the Rossiter-McLaughlin 
effect, using Eqs. J3J and jl| and assuming V ro tsin/* = 6.9 km s and <j> = 
20°, as estimated for KH 15D. 



dial velocity measurements that were obtained when the sys- 
tem flux was 90% or greater of its mean out-of-eclipse flux. 
This is a cautious procedure, in the sense that the systematic 
error is smaller for the more realistic case of a star with limb 
darkening and an occulting edge that is not perfectly sharp. 1 1 
We made two other minor adjustments to the data published 
by Johnson et al. (2004). First, we corrected the time stamps 
on some of the Keck measurements, which were in error by 
0.5 day because of a confusion between Julian and Modified 
Julian dates. Second, we added 0.25 km s -1 in quadrature to 
the uncertainties of the data labeled "b" in Table 1 of John- 
son et al. (2004). In those cases, the radial velocity scale was 

11 We also tried using all the data with / > 0, after enlarging the error 
bars to encompass the expected systematic errors. In that case, the procedure 
described in §|4]resulted in significantly poorer fits, but the optimized values 
of the parameters were similar to those reported in Table 1 . 
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TABLE 1 

Radial Velocity Measurements of 
KH 15D 



Julian Date 


Relative flux (/) 
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■i iO-U0^9 


2452653.9689 


0.6 


3.30 ±0.20 


2452678.9007 


0.7 


9.00 ±0.56 


2452679.9109 


0.2 


11.50 ±0.56 


2452707.5117 


1.0 


0.80 ±0.32 


2452947.0845 


0.9 


1.20 ±0.47 


2452948.0948 


0.9 


1.80 ±0.39 


2453008.7620 


1.0 


1.40 ±0.65 


2453009.8520 


1.0 


1.70 ±0.47 


2453014.7600 


0.9 


5.80 ±0.47 


2453014.8208 


0.9 


5.40 ±0.30 


2453016.0063 


0.4 


7.00 ±0.32 


2453044.8341 


1.0 


1.80 ±0.20 


2453045.8282 


1.0 


1.25 ±0.20 



NOTE. — From Hamilton et al. (2003) and 
Johnson et al. (2004), modified as described in 
§0 Positive velocity corresponds to motion 
away from the Sun. 



calibrated using only one reference star and the quoted error 
estimate was statistical, i.e., based only on the noise in the 
spectrum of KH 15D. In the other cases (labeled "a"), mul- 
tiple reference stars were available, and it was found that the 
statistical error was smaller than the variance in the results of 
using different choices for the reference star. Using the "a" re- 
sults we estimated that the systematic error in the calibration 
is w0.25 km s" 1 . For convenience, the complete and updated 
list of radial velocity measurements is given in Tabled Only 
the 12 measurements with / > 0.9 were used in the model- 
fitting procedure. 

4. THE MODEL OF KH 15D 
4.1. Definition of Parameters 

Our model consists of two stars (A and B) in a Keplerian 
orbit, and a straight-edged semi-infinite occulting screen rep- 
resenting the sky projection of a circumbinary disk. Star A 
is the star that is currently undergoing periodic eclipses, and 
for which radial velocity variations have been measured. The 
stars have masses denoted by Ma and Mb, and radii Ra and 
Rb- The orbit is specified by the period (P), eccentricity (e), 
argument of pericenter (lu), inclination with respect to the sky 
plane (/), heliocentric radial velocity of the center of mass (7), 
and a particular time of pericenter passage (T p ). 

A diagram of the coordinate system 12 is given in Fig. |4] 
The X-Y plane is the sky plane. The X-axis is chosen to be 
along the line of nodes, with the ascending node of star B at 
X > 0. The location and orientation of the occulting screen 
are specified by its point of intersection between the edge and 
the F-axis (Ye), and by the angle that the edge makes with the 
X-axis (9e)- In Paper I, we assumed that the screen moves 
uniformly and without rotation, i.e., Ye was taken to be a lin- 
ear function of time, and Be was taken to be a constant. In 
§|5]we consider these same assumptions, but we also consider 

12 Although our model is the same as that presented previously in Paper I, 
we have chosen a different coordinate system for convenience. 



more general cases that we argue are more realistic. 
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FIG. 4. — The coordinate system. The X-Y plane is the sky plane and 
the Z-axis points towards the observer. The orbital plane is inclined by an 
angle / with respect to the sky plane. The line of nodes is the X-axis, with 
the ascending node of star B occurring at X > 0. The Z > portions of the 
orbit are marked with solid lines, and the Z < portions are marked with 
dashed lines. The stars are shown at pericenter, and w is the argument of 
pericenter. The instantaneous location of the occulting edge is specified by 
Ye, its intersection with the 7-axis, and by 6e, the angle that the edge makes 
with the X-axis. A configuration with / < 90° , Ye < 0, and 6e < is shown. 

Once the orbital parameters and the screen's trajectory are 
specified, five "orbital contact times" can be defined (see 
Fig-H}. First contact (t\) occurs when the occulting edge is 
tangent to the sky-projected orbit of star B for the first time. 
Before first contact, no eclipses occur. Fig. |2] suggests that 
first contact occurred between 1959 and 1967. Second con- 
tact (?2) occurs when the occulting edge is first tangent to the 
sky-projected orbit of star A. Between t\ and t2, star B is pe- 
riodically occulted, but we receive constant light from star A. 
The total light exhibits "partial" or "diluted" eclipses, as were 
observed between the late 1960s and the early 1980s. After 
?2, both stars undergo periodic eclipses, and the light curve 
can appear quite complex. Third contact (tj) occurs when the 
edge crosses the projected center of mass of the binary. There 
are no abrupt changes in the photometric behavior. After ti, if 
the eccentricity is large, star B is seen for a rapidly shrinking 
fraction of the orbital period; it appears only briefly during 
the middle of each eclipse of star A. This is the era of "cen- 
tral re-brightenings," as observed in 1996 and 1997. After 
fourth contact (£4), the entire orbit of star B is hidden. As the 
screen advances further, the duration of the eclipses of star 
A increases until fifth contact (t^) when the entire system is 
covered. 

4.2. Interpretation of the Occulting Screen 

Of course, the real occulting edge is not semi-infinite. The 
dimensions of the screen's true sky projection depend on the 
three-dimensional structure of the circumbinary disk. Star B 
should eventually reappear, possibly even before t$, and the 
qualitative changes in total-light photometry that have been 
observed since 1960 may occur in reverse order as the trailing 
edge of the screen uncovers the binary. Those phenomena are 
not included in the model because the true three-dimensional 
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FIG. 5. — Definition of the 5 orbital contact times. In this illustration, the 
occulting screen has fefe) < 0, Be > 0, and Ye > 0. 



disk structure is not known and cannot be determined from 
existing observations. However, some remarks are in order 
regarding the relation between the idealized screen and a real 
circumbinary disk. 

The premise of the theory is that the circumbinary disk is 
inclined with respect to the binary plane, which causes the 
orbits of the disk's constituent particles to precess. Chiang & 
Murray-Clay (2004) showed that the disk should be warped, 
and that its inner and outer radii should be of the same order 
of magnitude, in order for the "ring" to precess as a unit. The 
edge of the screen may represent the ring's inner or outer edge 
in projection, or it may be a projection of the warp, in which 
case the edge does not necessarily represent a single orbital 
distance. These possibilities are depicted in Fig. [6] 

The correspondence between the precession of the disk and 
the motion of the screen is harder to visualize. For simplicity, 
imagine that the binary is surrounded by a circular ring of 
radius a r and inclination I r with respect to the binary's orbital 
plane. The ring represents, say, the inner radius of the disk, 
and its sky projection is the occulting feature. The ring's line 
of nodes (LON) precesses around the orbital plane with an 
angular frequency (3 < 0. 13 Now imagine viewing the stellar 
orbits edge-on. When the LON of the ring is pointed at the 
observer, the ring appears as a straight line with Y E = and 
9 E = I r - This is ?3, the moment of third contact. At other 
times, the sky projection of the ring is an ellipse. Considering 
the Z > half of the ellipse (i.e. the portion that is in front 
of the binary), the intersection point with the F-axis and the 
angle it makes with the X-axis are 



Y E (t) = a r 



e E (f)=wr l 



tan/ f sin/3(f-f3) 
yj 1 + sin 2 f3(t - 1 3 ) tan 2 I r 
tan/, cos f3{t-tx) [sin/3(f-f3)cos/tan/ r + sin/] 



cos 2 (3{t - 1 3 ) tan 2 (3{t - f 3 ) 



(5) 
4> 



For small /,., these expressions reduce to Y E (t) = 
a r tan/ r sin/3(f - t^) and 9 E (t) = 7, cos/3(f - t^). For times 

13 The negative value of /3 means that the LON regresses relative to the 
orbital angular momentum vector. 



close to f3, Y E (t) = /3a, tan/ r = constant, and 9 E (t) = /,• = con- 
stant, in accordance with the assumptions made in Paper I. 

However, the applicability of this approximation is not 
guaranteed. Over a 50-year time span, the higher-order terms 
in /3(t-tT,) may be appreciable. More importantly, the disk 
is probably elliptical, given that the stellar binary has e ~ 0.6 
(see Paper I, Hamilton et al. 2005, Herbst & Moran 2005, and 
also §[5]below). This alters both the amplitudes and phases of 
functions Y E (t) and E (t), depending on the eccentricity and 
argument of pericenter of the ring. These functions may also 
depend on the details of the warp. In addition, if a r is small 
enough, the curvature of the occulting edge in the sky plane 
may be significant. We found that there are too many possi- 
bilities for the data to distinguish among them. Instead, our 
approach was to begin by pursuing the simplest possibility, a 
straight-edged screen with a constant angle and speed, while 
also preparing to allow for angular rotation speeds (of order 
~/3/ r , according to the formula for a circular ring) and accel- 
erations (^/3 2 a r tan/ r ) if they appear to be needed. 

4.3. Expectations for Parameter Values 

Before proceeding to the model-fitting results, we review 
what was already known of the parameter values from other 
observations or theoretical arguments. Star A has been spec- 
trally classified as K6-7, with an effective temperature of ap- 
proximately 4000 K and luminosity 0.4 L & (Hamilton et al. 
2001, Agol et al. 2004), assuming a cluster distance of 760 pc 
(Sung et al. 1997). Its mass was estimated to be 0.6 Mq by 
Flaccomio et al. (1999) and 0.5-1.0 M by Hamilton et al. 
(2001), based on theoretical pre-main-sequence (PMS) evo- 
lutionary tracks. We find the range Ma = 0.6 ± 0. 1 Mq to be 
in reasonable agreement with the six different sets of PMS 
tracks compiled by Hillenbrand & White (2004). The corre- 
sponding radius estimate is Ra = 1.3 ±0.1 Rq. Unfortunately, 
we cannot estimate Ma and Ra independently using our model 
and the current data, because the limited radial velocity data 
place only weak bounds on the mass function of the binary, 
and because there is a degeneracy between the stellar radius 
and the projected orbital speed during occultations. 

No spectrum of star B has been obtained, but its luminosity 
can be estimated using simple considerations of the system's 
photometric history. The re-brightenings observed between 
1995 and 1998 (when star B was seen alone) had a maximum 
flux of 1.4 times the flux of the out-of-eclipse signal (when 
star A was seen alone), indicating that Lb /La ~ 1 .4. The pho- 
tographic photometry is not as accurate, but it also suggests 
Lb I La ^ 1. The flux of the maximum-light phase of the pre- 
1980 photometry (when both stars were seen) was approx- 
imately double the flux of the minimum-light phase (when 
star A was seen alone), and was approximately 2.3 times that 
of the modern maximum-light phase. The PMS tracks of 
D'Antona & Mazzitelli (1997) predict a mass-luminosity re- 
lation L oc M 19 , and a mass-radius relation R oc M° 21 , leading 
to the expectations M B /M A ~ 1.2 and Rb/Ra ~ 1-05. 

Johnson et al. (2004) analyzed the radial velocity data 
alone. Given the necessarily limited phase coverage of the 
observations, there were strong degeneracies between many 
of the orbital parameters. Among the most robust conclusions 
were that the orbital period is P = 48.4 days, the eccentric- 
ity is appreciable (e > 0.3), and the argument of pericenter is 
within about 20° of zero. The period determination is nearly 
identical with those based on periodograms of the photometry 
(by, e.g., Hamilton et al. 2005, Johnson & Winn 2004, Maf- 
fei et al. 2005), although we note that the periodogram-based 
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FIG. 6. — An illustration of the binary and the surrounding disk (above), and some possible viewing geometries (below). The true three-dimensional disk 
structure is not known. In particular the warp may have the opposite sign of the inclination gradient (dl r /da r ) as the warp depicted. In the lower panels, for 
clarity, only the cross-section of the disk in the sky plane is shown. 



method may be biased because the light curve is not strictly 
periodic. 

Herbst & Moran (2005) and Hamilton et al. (2005) made 
a theoretical prediction for the eccentricity, based on the the- 
ory of tidal interactions in close binaries. Tides will "pseudo- 
synchronize" the stellar spins and orbit, meaning that the stel- 
lar angular-rotation frequency, Q r , and the orbital angular fre- 
quency at pericenter, n p , will reach a theoretically determined 
equilibrium (Hut 1981): 



l + ip +fe 4 + ^6 



n p (l + 3e 2 +|e 4 )(l+e) 2 



^rb 
^rot 



(l-e 2 ) 3/2 
(1 + e) 2 



(7) 



The time scale for this process is poorly known but is es- 
timated to be a few million years, which is approximately 
the age of the system. It therefore seems likely that pseudo- 
synchronization is well underway, even if it has not yet been 
achieved. Hamilton et al. (2005) measured a rotation period 
P I0t = 9.6 days for star A, from periodic photometric variations 
observed out of eclipse and attributed to star spots. Given 
that P or b = 48.4 days, the pseudo-synchronization condition is 
e = 0.65±0.01.It seems likely that the rotation of star A has 
been slowed down, rather than sped up, given that its period is 



longer than usual for a low-mass star in NGC 2264. If pseudo- 
synchronization has not yet been achieved, then the equilib- 
rium period must be even longer, and the eccentricity must be 
correspondingly smaller to provide a slower angular velocity 
at pericenter. Hence, whether or not pseudo-synchronization 
has been achieved, the upper limit on the eccentricity is 0.66. 
This assumes Hut's theory to be applicable, i.e., the effects 
of dynamical tides, and of forces from the circumbinary disk, 
accreting material, or other bodies in the system, can be ne- 
glected. 

Some circumstantial evidence suggests that the orbital in- 
clination is nearly 90 degrees. First, the mere observation 
of eclipses suggests that the system is viewed nearly edge- 
on, unless the circumbinary disk is grossly misaligned with 
the orbital plane. Second, the inclination of star A's rotation 
axis (/*) seems to be near 90 degrees. The combination of the 
measured rotation period (P rot = 9.6 days), the projected rota- 
tion speed of the photosphere (V rot sin./* = 6.9 ± 0.3 km s -1 ), 
and the estimated stellar radius (Ra = 1.3±0.1 Rq) gives 
sin/* = 1 .01 ±0.09, or /* > 67°. Third, Hamilton et al. (2003) 
observed a double-peaked emission line of [O I], with the 
peaks at radial velocities ±20 km s" 1 of the systemic velocity. 
These are thought to arise from bipolar jets with space veloc- 
ities that are typically ±200 km s" 1 , suggesting that the jets 
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TABLE 2 

Observed Bounds on Eclipse Fractions 



Star 


Year 


Lower Bound 


Upper Bound 


g 


1958 


000 


050 


g 


1968 


282 


548 


g 


1970 


211 


634 


g 


1981 


275 


845 


B 


1996 


0.734 


0.957 


B 


1997 


0.920 


0.989 


A 


1997 


0.252 


0.360 


A 


1998 


0.270 


0.379 


A 


1999 


0.303 


0.569 


A 


2000 


0.330 


0.482 


B 


2000 


1.000 


1.000 


A 


2001 


0.326 


0.550 


A 


2002 


0.410 


0.432 


A 


2003 


0.444 


0.474 


A 


2004 


0.497 


0.539 


A 


2005 


0.530 


0.570 



(and by implication the stellar poles and the orbit) are inclined 
by > 80°. 

Finally, an expectation for 7 (the heliocentric radial veloc- 
ity of the center of mass) is available thanks to the radial- 
velocity survey of NGC 2264 by Soderblom et al. (1999). 
These authors found a median radial velocity of 20 ± 3 km s" 1 
among probable cluster members. Thus, although Johnson et 
al. (2004) found the radial velocity measurements to be com- 
patible with any value of 7 between 7 and 22 km s , we ex- 
pect the true value of 7 to be near the high end of that range, 
assuming KH 15D to be a cluster member. 

5. MODELS OF THE ORBIT 

We take a two-step approach to the quantitative modeling 
of KH 15D. The first step is a direct continuation of Pa- 
per I. Rather than fitting the photometry directly, we attempt 
to match some of its most important properties: the eclipse 
durations and the ingress and egress durations. For the mo- 
ment, we ignore all of the other information in the photome- 
try, such as the light observed during eclipses, and the detailed 
shapes of the ingress and egress light curves. The second step, 
described in § |6] is an attempt to fit the full photometric time 
series. 

To measure the eclipse durations, we plotted the phased 
light curve for each year in which fairly complete phase cov- 
erage was achieved, and determined the minimum and max- 
imum phases of the half-flux points of the eclipses. The re- 
sults are given in Table [2] They agree with those given in 
Paper I, and include new measurements based on the newly 
available photometry, as well as a few measurements of the 
eclipse fraction of star B based on the observations of central 
re-brightenings (or lack thereof). 

The ingress and egress durations are harder to measure be- 
cause fine time-sampling is needed. In Paper I we estimated 
the durations based on phased light curves from 2001-2003, 
but we see now that there is too much variation between 
eclipses for this method to be reliable. Instead, we used data 
from a single eclipse in 2002.1 that was observed with espe- 
cially high cadence. We considered only the steepest (most 
rapidly varying) portion of the light curve between flux levels 
33% and 67% of the maximum, and fitted a model light curve 
to these data consisting of a limb-darkened star being crossed 



by a knife-edge occultor. The best-fitting ingress and egress 
durations (the times taken to cross the diameter of the star) 
were 2.55 ± 0.08 days and 2.96 ± 0.05 days, respectively. 

We assume Ma = 0.6 M Q and Ra= 1-3 Rq as explained in 
§ 14.31 The mass of star B is taken to be a free parameter. There 
are 6 free parameters describing the orbit: {P,e,I,ui,'y,T p }. 
Additional free parameters are needed to describe the occult- 
ing screen and its trajectory. We begin with the case of a 
screen that moves uniformly with a fixed angle, in which case 
3 parameters are needed. One obvious parameterization, for 
example, is {Ye(T p ),Ye,0e}- We use the alternative parameter 
set {f4,?5, 9e}, which provides a more direct connection to the 
observations. 

The figure-of-merit function is 

222 
X =X V + Xd 

where vo,i are the observed radial velocities (of which there 
are N v = 12), <7 V)I are the corresponding 1 a uncertainties, and 
Vc,i are the radial velocities as calculated according to the 
model. The single measurement of the ingress duration (dc,\) 
and the egress duration (dc 2) are described with a similar no- 
tation. We used an AMOEBA algorithm (Press et al. 1992, 
p. 408) to minimize x 2 , while simultaneously requiring the 
calculated eclipse durations to obey the 16 constraints 14 given 
in Table 13 We also enforced the proper phase alignment be- 
tween the radial velocity data and the photometry by requiring 
the model system to experience mid-eclipse within 5 days of 
JD 2,452,352.5, a particular mid-eclipse observed in 2002. 

An excellent fit is achieved, with % 2 = 3 and 4 degrees of 
freedom. Thus the data are simultaneously consistent with 
the measurements of radial velocity, ingress and egress dura- 
tion, and eclipse duration. FigureQshows the fit to the radial- 
velocity and eclipse-fraction data. It also shows a diagram of 
the orbit and the screen at each of the 5 orbital contact times. 
The best-fitting parameter values are given in Table [3] under 
the heading "Model 1." 

Many of the best-fitti ng p arameter values conform with the 
expectations given in § 14.31 the period is P = 48.4 days, the 
eccentricity is e = 0.55, the inclination is 83°, and the argu- 
ment of pericenter is 6°. However, there are two exceptions. 
First, the heliocentric radial velocity of the center of mass (14 
km s" 1 ) is somewhat low; KH 15D would be a ^2 a out- 
lier in the radial velocity distribution of the cluster. A sec- 
ond and much more serious problem is that the mass ratio is 
Mb I 'Ma = 0.76, which is less than unity, even though it seems 
certain that Lb /La > 1 . 

The latter problem can be traced to the assumption that the 
occulting screen moves uniformly in a constant direction. Un- 
der this assumption, ts —ti is proportional to the projected size 
of star A's orbit, which in turn is proportional to M~ A X . Like- 
wise, t\-t\ oc Mg l . The photometry provides the constraints 
1959 < h < 1967, 1981 < t 2 < 1996, 1997 < 1 4 < 2000, and 
ts > 2005, and an extrapolation of the rapidly rising eclipse 

14 The eclipse duration measurements have the character of strict bounds, 
rather than central values and 1 a uncertainties, because the measurements 
are limited by time-sampling of the events rather than statistical errors in the 
flux measurements. 
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FIG. 7. — Model 1: Mb is a free parameter and 8e = 0. Left. — The binary orbit (solid lines) and the position of the occulting edge at the 5 orbital contact times 
(dashed lines). Middle. — The eclipse fractions of stars A and B as a function of time, as observed (error bars) and calculated (solid lines). Right. — The radial 
velocity variations of star A as a function of orbital phase, as observed (points with error bars) and calculated (solid line). 



fraction suggests t$ < 2010. Together, these inequalities give 
M B t 5 - 1 2 

— = < 0.97 (assuming constant screen velocity). 

M A f 4 -fi 

(9) 

If Mb /Ma is held fixed at 1.2, the value expected from a con- 
sideration of the luminosity ratio, then it is impossible to fit 
the eclipse data. For example, in the model, star B is covered 
6 years too early, or the rise in the eclipse fraction of star A 
over the last decade is m uch slower than has been observed. 

As explained in § 14.21 one might expect the real projection 
of a circumbinary disk to exhibit deviations from a constant 
velocity. Once the screen is allowed to rotate or accelerate 
during its passage in front of the binary orbit, a much wider 
range of mass ratios is allowed. While this saves the model 
from a seemingly unphysical prediction, it becomes impossi- 
ble to determine the stellar mass ratio independently. Instead, 
in what follows, we hold the mass ratio fixed at Mb /Ma = 1 .2. 
Because the leading-order correction in the precession angu- 
lar velocity (3 is a rotation, we concentrate on the case of uni- 
form motion plus a uniform rotation (with no acceleration). 
After dropping one of the original free parameters affecting 
the radial velocities (Mb), and adding a new one (8e) affecting 
only the eclipses, the resulting fit is also good, with x 2 = 6 and 
5 degrees of freedom. The fit to the data is shown in Fig. [8] 
and the best-fitting parameter values are in Table[3]under the 
heading "Model 2." 

We conclude that it is possible to fit the data with a phys- 
ically plausible mass ratio. All of the parameters agree with 
the expectations given in § 14.31 including the heliocentric ra- 
dial velocity (7 = 18.4 km s" 1 ). The best-fitting value for the 
rotation rate is 8e = 6 x 10~ 3 rad yr -1 , which is compatible 
with reasonable guesses for the circumbinary disk properties. 
It implies a precession rate of order (3 ~ 10~ 2 rad yr" 1 , as 
would be experienced by a circular ring of radius ~3 AU. 
However, a perfectly circular ring is not realistic on physical 
grounds, and is not consistent with the data beyond this order- 
of-magnitude comparison; the rotation rate in that case would 
be exactly zero at t-$, in contradiction to the assumptions of the 
model. We also found that a uniform acceleration with no ro- 
tation provides a reasonable fit, and obviously a combination 
of rotation and acceleration is possible. Given the uncertainty 



in both the stellar mass ratio and the geometry of the disk, we 
cannot use the model alone to decide which of these cases is 
closer to the truth. 

6. MODELS OF THE OCCULTATIONS 

The models described in the previous section ignore some 
potentially interesting and information-bearing aspects of the 
photometry, such as the phase of each individual eclipse, the 
shape of the ingress and egress light curves, and the progres- 
sive deepening of the eclipses. In this section, we elaborate 
upon our model in several ways and attempt to fit all of the 
photometry. 

In Paper I, we generated model light curves by treating the 
occulting feature as a completely opaque and semi-infinite 
knife-edge, and the stars as limb-darkened photospheres. Al- 
though these model light curves successfully reproduce the 
basic patterns in the photometry, there are at least three ways 
in which they fail to match the data within the measurement 
uncertainties: 

1 . The model does not reproduce the observed ~0. 1 mag 
variations outside of eclipses. As noted previously, 
Hamilton et al. (2005) provided evidence for quasiperi- 
odic variability outside of eclipses, and interpreted the 
9.6 day period as the stellar rotation period of star A. 
Presumably, these episodic variations are caused by 
cool star spots that last only a few months, as is com- 
mon for young low-mass stars (see, e.g., Herbst et al. 
1994). 

2. The steepest (central) phase of each occultation light 
curve is consistent with a knife-edge crossing a photo- 
sphere. However, during the first and last third of each 
event, the flux variation is slower than such a model 
would predict, as noted by Agol et al. (2004) and in 
Paper I. 

3. Light is observed even during mid-eclipse, when (ac- 
cording to the model) the photospheres of both stars 
are completely hidden. The mid-eclipse flux has fallen 
from about 9% to 2% of the uneclipsed flux over the 
last decade (see Fig.|2jl. Indeed, the mid-eclipse light is 
not constant even within a single eclipse. Although the 
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FIG. 8. — Model 2: Mb/Ma = 1.15 and Oe is a free parameter. The plotting conventions are the same as in Fig. 171 



re-brightening events have diminished dramatically in 
intensity since 1998, there are still faint echoes of this 
phenomenon in the more recent data. 

We deal with the first point by treating the episodic varia- 
tions of star A as an intrinsic source of noise. We do not ex- 
pect the light curve to be predictable with better than 0.1 mag 
accuracy. As for the second and third points, we have found 
that a good fit is achieved with a simple model in which each 
star is enveloped by a faint halo. Fig. [9] illustrates how this 
model accounts for the the light observed during eclipses, 
the intra-eclipse variations, and the shape of the ingress and 
egress light curves. The dashed lines mark 6 special photo- 
metric phases, labelled a—f, for which the configurations of 
the photospheres, the halos, and the occulting screen are il- 
lustrated in the smaller panels beneath the light curve. These 
phases are: (a) Out of eclipse. Both the photosphere and the 
halo of star A are seen, (b) Start of ingress. The halo of star 
A begins to be occulted, (c) Midpoint of ingress. The pho- 
tosphere is half-covered and the total flux declines steeply, 
(d) The photosphere is completely occulted, and the total flux 
resumes its shallower decline, (e) The halo of star A is al- 
most completely covered, and the total flux reaches a local 
minimum, (f) Star B reaches its closest approach to the oc- 
culting edge, and its halo is maximally exposed. The total 
flux reaches a local maximum. (At this phase, before U, the 
photosphere of star B was exposed and manifested as the re- 
brightening event.) Although it is not illustrated in Fig. |9] the 
gradual darkening of the mid-eclipse phase can also be ac- 
counted for with this model: as the occulting edge advances, 
the fraction of the halo of star B that is exposed at mid-eclipse 
[panel (f)] gradually shrinks. 

To develop a quantitative model, we first demonstrate that 
the occultations of star A are repeatable events: all of the 
observed ingress and egress events are nearly equivalent, after 
correcting for the differences in orbital velocity of the star at 
the contact points with the occulting edge. Using our model, 
we calculate the position of star A and of the occulting edge 
as a function of time. Then, instead of plotting the flux vs. 
time, we plot the flux vs. the position of the occulting edge 
relative to star A. For this purpose it is convenient to use a 
rotated coordinate system (x,y), in which the occulting edge 
is the line x = xg(f), the perpendicular direction is the y-axis, 



and distances are measured in units of Ra- The sky position 
of star A is (xa,va). The top panel of Fig.llOlshows the result, 
summarizing 10 years of photometry as star A goes back and 
forth beneath the occulting edge. The measured flux is indeed 
a coherent function of Ax = x E -xa, with a scatter of ~ 10%. 

This allows us to model the occultations as a knife edge 
scanning over a time-invariant surface brightness distribution 
associated with the star, S(x-XA,y-yA). Since the relative 
velocity of the edge and the star is always closely aligned with 
the x-axis, the only relevant quantity is the one-dimensional 
brightness distribution 



B(x)-. 



dy S(x,y). 



(10) 



When the occulting edge is at xe and star A is at (xa^a), the 
measured flux is the cumulative 1-d brightness distribution 



C = 



dx B(x—xa) = 



dx .B(x), 



(11) 



which is a function of Ax, as observed. We note that B rep- 
resents all of the flux that moves with star A on the sky. It 
includes light from the star, any luminous material that orbits 
with the star, and also any scattering halo (which we will later 
argue does exist). The halo need not be physically associated 
with the star, just as the lunar halo that is sometimes observed 
on cold nights is not physically associated with the Moon. 

Next, we devise a fitting function for B that is consistent 
with the top panel of Fig. EH At Ax = 0, the occulting edge 
intersects the center of star A, and approximately one-half of 
the total flux is exposed. As the screen advances over the star 
from Ax = to 1 , the received flux falls rapidly, in a manner 
consistent with a knife edge covering a stellar photosphere. 
The abrupt slope change at Ax = 1 is suggestive of the contact 
between a sharp occulting edge and the rim of the stellar pho- 
tosphere. For Ax > 1, the flux falls off much more gradually, 
as if the star is trailed by a spatially extended halo comprising 
a few per cent of the total flux. Interestingly, the light curve 
is not symmetric about Ax = 0. Once the star is more than 
half exposed (Ax < 0), the flux rises more gradually than one 
would expect for a symmetric halo. The side of the star near- 
est to the occulting edge appears more smeared out than the 
other side. 
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FIG. 9. — Illustration of the model for the occultation light curve. Top. — Phased light curve from the 2001-2002 season, with 6 particular phases marked with 
vertical dashed lines. Bottom panels. — Corresponding configurations of the stars, halos, and the occulting edge. These are cartoons only, and do not represent 
optimized model parameters. In particular, the best-fitting model halos are asymmetric unlike the circular halos drawn here. 



After some experimentation, we found good results for the 
fitting function 

fBiexp[(x+l)/ei], x<-\ 
B(x) = < B x +B±(x), -1 < x < 1 (12) 

[B 2 exp[-(x-l)/6] , x>l. 

The ingredients of this function are the brightness distribu- 
tion of the stellar photosphere (fi*), and an exponential fall- 
off on each side of the star. The amplitude (B\ or B2) and 
scale length (£1 or £2) of the exponential function are allowed 
to be different on each side, to match the observed asymme- 
try in C, and the function is discontinuous at x = 1, to match 
the observed slope discontinuity in C. In practice, rather than 
quoting B\ and B2, it is more useful to describe the bright- 
ness distribution by the total flux (L) and the fraction of the 
total flux carried by each side of the exponential halo (ei and 
62). Although the photosphere may be modeled as uniform 
to within the accuracy of the data, for completeness we use 
a linear limb-darkening law, with a coefficient of 0.65. 15 We 
emphasize that Eq.^]is merely a fitting formula. In fact, we 

15 For a star with r eff 4000 K and g = 10 4 cms -, a linear limb-darkening 
coefficient of u fa 0.65 is predicted by Claret ( 1 998) and Claret (2000) for the 
Cousins / band. Van Hamme (1993) predicts a similar value, 0.6. The value 
of 0.3 used by Herbst et al. (2002) and in Paper I was chosen mistakenly. In 



will argue later that the halo is not a physical object, but rather 
represents forward- or back-scattering by particles in the cir- 
cumbinary disk. In this section, we limit ourselves to demon- 
strating the success of the phenomenological "halo model," 
deferring further discussion of the physical interpretation un- 
til §□ 

Next, we proceed to the optimization of parameters. As be- 
fore, we assume M A = 0.6 M Q , R A = 1.3 R & , and M B /M A = 
1.2. The occultor is an opaque, semi-infinite straight edge 
whose intersection point Ygit) and rotation angle £?£(f) are 
both linear functions of time. Both star A and star B are 
modeled with the function given in Eq. (I12> . and we add a 
time-invariant term Lq, representing any large-scale (unoc- 
culted) components of the surface-brightness distribution. All 
together, the model has N p = 17 free parameters: 7 specifying 
the surface-brightness distribution, {La,Lb,Lq,ei, £2,£,i,£,i}', 
6 specifying the orbit, {P,e, I,w,T p , 7}, and 4 describing the 
occulting screen, {f^fs,^^),^}. 

It was desirable to speed up the computations by time- 
averaging the photometry. Whenever more than one measure- 
ment was made within a given 6 hour time period, we com- 

any case, the choice of the limb-darkening coefficient has a very minor effect 
on the model light curve. 
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FIG. 10. — 7b/>. — The relative flux of KH 15D as a function of the position 
of the occulting edge. The quantity A.v is the sky-projected distance between 
the edge and the center of star A, in units of the stellar radius. For a perfectly 
sharp and straight edge, this plot can be interpreted as the cumulative 1-d 
brightness distribution of the star and its apparent halo. The dots show the 
data, and the solid line is the best-fitting model. Bottom. — The 1-d brightness 
distribution of the best-fitting model. 



TABLE 3 
KH 15D Model Parameters 



Parameter Best-fitting Value 





Model 1 


Model 2 


Model 3 


M A [M Q ] 


0.6 


0.6 


0.6±0.1 


M B /M A 


0.76 


1.2 


1.2±0.1 


L A [Lq] 






0.413 ±0.006 


Lb /La 






1.36 ±0.10 


Lo/L A 






0.029 ± 0.005 


Ra [Rq] 


1.3 


1.3 


1.30 ±0.07 


Rb/Ra 


1.05 


1.05 


1.05 ±0.03 


Si 






2.0 ±0.3 








3.3 ±0.4 








0.25 ±0.05 


C2 






0.049 ± 0.003 


i [Udysj 


48 4 


48 4 


to.joi in u.wj 


e 


0.55 


0.58 


0.574 ±0.017 


I [deg] 


83 


83 


92.5 ±2.5 


w [deg] 


6 


14 


13±2 


T p [JD] - 2,452,350 


1.1 


0.7 


1.9±0.7 


7 [km s -1 ] 


14.1 


18.4 


18.6=1= 1.3 


t\ 


1965 


1963 


1962 ±3 


h 


1985.3 


1985.9 


1985.9 ±0.7 


h 


1991.1 


1992.5 


1992.7±0.5 


n 


1998.8 


1997.0 


1997.3 ±0.5 


H 


2010.9 


2008.0 


2008.0 ±0.2 


9e(.U) [deg] 


-20 


-16 


-32 ±7 


6 E [radyi- 1 ] 





0.0061 


0.0087 ±0.0015 



NOTE. — Models 1 and 2 were fitted to the eclipse du- 
ration, ingress/egress duration, and radial velocity data. In 
Model 1, M B is a free parameter and 9 E = 0. In Model 2, 
Mb I Ma =1.2 and 8e is a free parameter. Model 3 was fit- 
ted to the photometric and radial velocity data. The uncer- 
tainties were determined with the Monte Carlo procedure 
described in the text. 



puted the flux-weighted mean of the results. The resulting 
binned data set has Nj = 1044 entries. The fitting statistic is 




which is analogous to Eq. (|8j, with the fluxes / replacing the 
ingress or egress durations d. The uncertainty oy ,- was taken 
to be the quadrature sum of the measurement uncertainty and 
10%, where the latter is intended to account for intrinsic vari- 
ations. The multiplier A controls the relative weighting of the 
flux data and radial velocity data. If the measurement un- 
certainties were known perfectly, and the model were correct, 
then A = 1 would be the appropriate choice and the best-fitting 
solution would have \ 2 w 1040 = Nf+N v -N p . 

In this case, A = 1 is probably not the best choice. We expect 
the model to be statistically consistent with the radial veloc- 
ity data, because the only relevant aspect of the model — the 
assumption of a fixed Keplerian orbit — seems reasonable. In 
contrast, the photometric model involves idealized assump- 
tions about the occultor and the surface-brightness distribu- 
tion, and the estimates of ay are crude. Our approach was to 
force the model to be statistically consistent with the radial ve- 
locity data, by increasing A until the best-fitting solution had 
Xv ~ 12 = /Vy . Of course, we cannot be certain that the radial- 
velocity data are free from additional sources of measurement 
error and unmodeled systematic effects, such as precession of 
the orbit due to perturbations from the circumbinary disk and 
additional bodies. Nevertheless, in the rest of this section, we 
describe the optimized model that results from this procedure. 



With A = 50, the best-fitting model has \v = 13 and xj = 
1574. The orbit, and the fit to the radial velocity data, are 
shown in Fig. The fit to the photometry is shown in 
Fig. Each panel shows the phased light curve for a par- 
ticular time interval. The time intervals chosen for this plot 
are the intervals that are marked by dashed lines in Fig. |2] 
The final column of Table [3] under the heading "Model 3," 
gives the best-fitting parameters. 

The table includes estimated uncertainties in the fitted 
quantities, which should be treated with caution. They were 
calculated using a Monte Carlo algorithm, as follows. We 
optimized the parameters on each of 10 4 synthetic data sets 
("realizations"). Each realization consisted of Nf = 1044 flux 
measurements and /Vy = 12 radial velocity measurements (the 
same numbers as the real data set). The flux measurements 
were randomly selected from the real data set, with repeti- 
tions allowed. The intent is to estimate the probability distri- 
bution of the measurements using the measured data values 
themselves (see, e.g., Press et al. 1992, p. 689). The num- 
ber of radial velocity measurements is too small for drawing- 
with-replacement to be effective. Instead, we added Gaus- 
sian random numbers to the measured radial velocities, with a 
mean of zero and a standard deviation given by the quoted a v 
values. To account for the uncertainty in the stellar masses, 
we assigned Ma/Mq by picking a random number from a 
Gaussian distribution with mean 0.6 and standard deviation 
0.1. Likewise, for the mass ratio Mb /Ma, we used a Gaus- 
sian distribution with a mean of 1.2 and a standard deviation 
of 0.1. The stellar radii were fixed according to the relation 
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R oc M . For each parameter, we found the median and the 
approximate 68% confidence limits in the resulting distribu- 
tion of best-fitting values. 

There are many caveats. To begin with, the model is not 
in formal statistical agreement with the data (some discrepan- 
cies are highlighted below). The random flux variations of the 
stars are not independent; they are correlated on time scales 
of days. The a priori probability distributions for the stellar 
masses are probably not Gaussian. The systematic error re- 
lated to our assumption about the trajectory of the occulting 
edge has not been taken into account. One can imagine more 
complex analyses that attempt to overcome these problems, 
but we do not believe that such methods are justified. 

The model shares the successes of the models presented in 
§ 13 It reproduces the transitions from no eclipses to partial 
eclipses to total eclipses, the existence of the central reversals 
from at least 1995 until 1999, and the growing duty cycle of 
the modern eclipses. The more detailed photometric model 
also accounts for the relative phases of the partial and total 
eclipses, the gradual deepening of the eclipses and the slope 
discontinuities in the ingress and egress light curves. In addi- 
tion to fitting the data reasonably wel l, the model satisfies all 
of the criteria that were discussed in § 14.31 The stellar masses, 
radii, and luminosities are reasonable, by fiat. The orbital ec- 
centricity is smaller than the theoretical upper bound of 0.66, 
and the orbital inclination is close to edge-on. The heliocen- 
tric radial velocity of the center of mass is near the median 
that has been observed among cluster members. The screen 
trajectory is consistent with order-of-magnitude estimates for 
the expected motion of the edge of a precessing circumbinary 
ring with a radius of ^3 AU. 

The best-fitting brightness distribution for star A is shown 
in Fig.^3 The upper panel shows the model cumulative 1- 
d brightness distribution overlaid on the data points, while 
the lower panel shows the model brightness distribution it- 
self. The asymmetry of the halo is evident. The exponential 
scale length on both sides is a few stellar radii, but the frac- 
tional flux of the halo is much greater on the side closest to 
the occulting edge (~25%) than on the other side (~5%). 

As noted before, the model cannot reproduce the apparently 
random ~10% flux variations of star A that occur on time 
scales of a few days. Nor can it reproduce the seasonal fluc- 
tuations of a few per cent in the mean flux of star A. There 
are some other significant discrepancies between the data and 
the model. In the 1967-1972 time period, there is a ~0.1 mag 
offset between the calculated and observed light curves, and 
in addition, the calculated eclipse depth is ~0.15 mag too 
deep. These problems are at least partly the fault of the un- 
certainty in the photographic magnitude scale. The offset is 
within the 0.14 mag uncertainty in the zero point, and the cal- 
culated depth becomes too sha llow if we use the Im system 
instead of the Ij system (see § 12. 1> . Another discrepancy is 
the phase shift of ^0.02 between some of the calculated and 
observed central re-brightenings. This is most easily seen in 
the 2000-2001 data. This offset disappears for an orbital pe- 
riod of 48.36 days, but that value of P is disfavored by the 
radial velocity data and by the relative phases of the pre- 1980 
eclipses. Finally, there is one data point in 1997-1998 that 
is very poorly described: at / = 16.5, it would seem that the 
photosphere of star B was partly exposed, but in the model, it 
was behind the screen. 



7. SUMMARY AND DISCUSSION 

Based on the preceding results, we draw two conclusions. 
First, the newly-available data corroborates the theory that 
KH 15D is an eccentric binary system being occulted by the 
edge of a precessing circumbinary disk. As shown in §|5] the 
model presented in Paper I succeeds as a quantitative descrip- 
tion of KH 15D, even when confronted with far more data 
than were available when the model was invented. The best- 
fitting parameter values are realistic and conform to theoreti- 
cal expectations. A problem that was overlooked in Paper I — 
the seemingly unphysical mass ratio of the stars — is fixed with 
a minor and physically-motivated elaboration to the model, 
namely, the allowance for a more realistic trajectory for the 
occulting edge across the line of sight. Second, the occulta- 
tion light curves are well-fitted by a model in which each star 
is surrounded by a more extended halo. This model succinctly 
accounts for both the photometric variations observed during 
individual events, and the gradual deepening of the eclipses. 
The halo around each star is not symmetric about the center 
of the star, but in each direction it has a typical scale length of 
a few stellar radii, and it is brighter in the direction facing the 
occulting edge. Only the one-dimensional profile of the halo 
has been measured. 

What are these halos? The three broad categories of pos- 
sible physical interpretations of each halo are (1) luminous 
material, such as accretion columns or a hot corona, (2) par- 
tially extinguished starlight, (3) scattered starlight. The evi- 
dence does not permit an unambiguous interpretation, but it 
does strongly suggest that the halo light has a scattered com- 
ponent. Relative to the uneclipsed light, the mid-eclipse light 
has a larger fractional polarization (Agol et al. 2004) and is 
slightly bluer in color (Hamilton et al. 2001, 2005; Herbst et 
al. 2002; Knacke, Fajarado-Acosta & Tokunaga 2004; Kusak- 
abe et al. 2005). Enhanced polarization and blueness are hall- 
marks of scattering by small particles. In contrast, partially 
extinguished light would be redder than the light source (or 
the same color, if the particles causing the extinction were 
large enough). Likewise, the observation of only a small color 
change argues against self-luminous material in an accretion 
flow or a hot corona, which would not be expected to have 
nearly the same effective temperature as the photosphere. 

Just as the brightness distribution of the halo was deter- 
mined from /-band photometry in § [6] the color distribution 
of the halo can be measured via multi-band photometric mon- 
itoring. A first attempt at this is shown in Fig. [O] which is 
analagous to Fig. ^| but uses the multi-band photometry of 
Hamilton et al. (2005) and Kusakabe et al. (2005) rather than 
monochromatic photometry. Higher-precision monitoring in 
multiple bands throughout an ingress or egress will help to 
resolve the color structure of the halo. 

As for the location of the scattering material, the two basic 
possibilities are the immediate vicinity of each star, and the 
circumbinary disk. Examples of the first scenario are scatter- 
ing from a corona of hot electrons, or from infalling material. 
A difficulty with this scenario is the near-repeatability of the 
occultation light curve. A circumstellar halo would need to 
maintain a fixed orientation in space with respect to the oc- 
culting feature, despite the rotation of the star and its tidal 
interaction with its binary companion. Disk scattering, on the 
other hand, would need to be dominated by forward-scattering 
or back-scattering (or some combination of the two), in or- 
der to give the appearance of a localized halo around the star. 
The more general illumination of the disk would not depend 
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strongly on the positions of the stars, because the disk ra- 
dius of ^3 AU is much larger than the binary semimajor axis 
of ^0.25 AU. It is possible, however, that this is a false di- 
chotomy. The disk may extend inward continuously to much 
smaller distances, and indeed the strong Ha emission, indica- 
tive of ongoing accretion, can be taken as evidence for this. 

One way in which these possibilities might be distinguished 
is through the difference in the typical scattering angle. In the 
case of circumstellar scattering, the received light is scattered 
over a wide range of angles, as opposed to the forward- or 
back-scattering from the disk. The chromatic and polariza- 
tion signals associated with scattering depend on the typical 
scattering angle, as well as the size and composition of the 
particles. Hence, measurements of the polarized light curve, 
combined with models of the expected polarization signal un- 
der various circumstances, may help to determine the size, 
composition and spatial distribution of the scattering parti- 
cles. It would also be possible to determine the velocity of 
the scattering material relative to the stars, through observa- 
tions of Doppler shifts in the stellar absorption features of the 
scattered light. 

Another interesting feature of the model is the extreme 
sharpness of the occulting edge. The best evidence for a 
sharp edge is the slope discontinuity in the light curve that 
was pointed out in §[6] This "kink" is most easily understood 
as the point where the sharp occulting edge contacts the rim 
of the photosphere of the star. If the edge were not sharp, 
then the time of photospheric contact would be smeared out. 
The presence of the kink requires that the optical depth of the 
edge increases from nearly zero to >3.5 over a distance that 
is smaller than a stellar radius. It can also be argued that the 
edge is sharp on the scale of a stellar radius based on the dra- 
matic changes in the Ha emission and absorption profile that 
Hamilton et al. (2003) observed during an occultation. The 
implication is that the occulting edge spatially resolves the 
Ha-emitting region, which is thought to arise within a few 
stellar radii (i.e. within the magnetosphere). 

An independent way to assess the sharpness of the edge 
would be to measure the Rossiter-McLaughlin effect. In §[5J 
we predicted the amplitude of this effect for a perfectly sharp 
and straight edge. In the other extreme, in which the scale 



length of the edge is much larger than the stellar radius, there 
is no Rossiter effect. Thus, a sequence of optical spectra with 
a high signal-to-noise ratio throughout an occultation could be 
used to determine whether or not the occulting edge resolves 
the stellar photosphere. If the occulting edge is truly sharp 
on the scale of the stellar radius, then it can be used to map 
out the environment of star A on scales unobservable by any 
other means, just as lunar occultations were used to identify 
radio and X-ray sources before the necessary angular resolu- 
tion was available. Spectroscopic monitoring would resolve 
the velocity structure of the surrounding material. The pro- 
jected orbital velocity of star A relative to the edge is approx- 
imately one stellar radius per day. Thus, a time sampling of 1 
hour corresponds to a spatial resolution of ~0.05 stellar radii, 
or an angular resolution of ~0.5 ^as at a distance of 760 pc. 
However, the scattering geometry will need to be understood 
first, to allow the information in the halos to be decoded. 

An interesting theoretical question is why the edge is sharp. 
With a radius of a r = 3 AU and a vertical scale height h < Rq, 
the circumbinary disk would have h/a r < 10~ 3 . As pointed 
out by Chiang & Murray-Clay (2004), this is at least an or- 
der of magnitude smaller than the values generally assumed 
for T Tauri stars, which are based on the consideration of hy- 
drostatic equilibrium between gas pressure in the disk and the 
vertical gravitational field provided by the star. One possi- 
ble resolution is that there is little or no gas over the range 
of orbital distances that form the occulting feature. Another 
possibility is that a r is smaller than the dynamical arguments 
suggest. A third and most intriguing possibility is that the 
obscuring particles in the disk have settled to the midplane 
and formed a much sharper layer than the distribution of gas. 
The dust-settling process has long been investigated as a pos- 
sible mechanism for generating planetesimals, by enhancing 
the surface density of solid particles until a gravitational insta- 
bility ensues (Safronov 1969, Goldreich & Ward 1973, Wei- 
denschilling 1977). 

Finally, we end this discussion on a note of urgency. A 
problem with all of the observations proposed above is that 
the edge is continuing to advance and will eventually cover 
star A at all orbital phases. According to the model of §[6] the 
photosphere of star A will set for the last time near the begin- 
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FIG. 12. — Model 3: Calculated and observed light curves. Each panel shows the phased light curve from a particular time range, as observed (solid symbols) 
and calculated (dots). The boundaries of the time ranges are marked with dashed lines in Fig. [2] 
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ning of 2008, although that prediction hinges on our extrap- 
olation of the trajectory of the occulting screen. The halo of 
star A will be visible for several years after the photosphere is 
hidden. We encourage continued monitoring and will gladly 
provide predictions of future occultations to interested parties. 



*A) / *A 



FIG. 13. — The color index of KH 15D, as a function of the distance 
between the occulting edge and the center of star A. The median color out- 
of-eclipse is indicated with a solid line. To make this plot, the time stamps 
on the V — Ij photometry of Hamilton et al. (2005) and the J—H photometry 
of Kusakabe et al. (2005) were converted into (xe ~xa), using the model 
presented in g |S] The halo is seen to be bluer than the photosphere at both 
optical and near-infrared wavelengths. 
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